arXiv:astro-ph/0112495vl 20 Dec 2001 


Mon. Not. R. Astron. Soc. 000, 000-000 (0000) 


Printed 5 February 2008 (MN style file vl.4) 


Time evolution of the linear perturbations of a rotating 
Newtonian polytrope 

D. I. Jones^, N. Andersson^ and N. Stergioulas^ 

^ Faculty of Mathematical Studies, University of Southampton, Highfield, Southampton, S017 IBJ, United Kingdom 
^ Department of Physics, Aristotle University of Thessoloniki, Thessaloniki 54006, Greece 


5 February 2008 


ABSTRACT 

We present the results of numerical time evolutions of the linearised perturbations 
of rapidly and rigidly rotating Newtonian polytropes while making the Cowling ap¬ 
proximation. The evolution code runs stably for hundreds of stellar rotations, allowing 
us to compare our results with previously published eigenmode calculations, for in¬ 
stance the f-mode calculations of Ipser & Lindblom, and the r-mode calculations of 
Karino et al. The mode frequencies were found to be in agreement within the expected 
accuracy. We have also examined the inertial modes recently computed by Lockitch 
& Friedman, and we were able to extend their slow-rotation results into the rapid 
rotation regime. In the longer term, this code will provide a platform for studying a 
number of poorly understood problems in stellar oscillation theory, such as the effect 
of differential rotation and gravitational radiation reaction on normal mode oscilla¬ 
tions and, with suitable modifications, mode-mode coupling in the mildly non-linear 
regime. 

Key -words: stars: neutron - stars: rotation 


1 INTRODUCTION 

The study of the normal modes of stars has long been recog¬ 
nised as a subject of considerable importance. The electro¬ 
magnetic signatures of such oscillations have been used for 
decades to probe the structure of main sequence stars (Cox 
1980). More recently the modes of compact stars (neutron 
stars and white dwarfs) have received attention from the¬ 
orists, again with a hope of converting observations into 
constraints on stellar structure. Such modes may have elec¬ 
tromagnetic signatures (van Horn 1980), excited perhaps 
during accretion or in a thermonuclear burst (Van der Klis 
2000). Equally, an oscillating compact star is a source of 
gravitational radiation, and may even be unstable to such 
an energy loss via the so-called Chandrasekhar-Friedman- 
Schutz mechanism (hereafter CFS; see Friedman & Schutz 
1978a,b). Particularly in this last regard, it is rapidly ro¬ 
tating compact stars that are of most interest, as only in 
the case of rapid rotation might gravitational radiation re¬ 
action overcome the dissipative effects of viscosity to amplify 
a mode. 

Very few previous studies have examined modes of 
rapidly rotating compact stars in Newtonian theory (see 
Stergioulas 1998, Andersson & Kokkotas 2001 for reviews 
of relativistic computations). The f-modes were calculated 


by Ipser & Lindblom (1990), the r-modes by Karino et al. 
(2000), and the inertial modes by Lockitch & Friedman 
(1999). All of these calculations were carried out by means 
of an eigenmode analysis, i.e. by assuming an oscillatory 
constant-frequency solution. In contrast, in this paper we 
will perform time evolutions of the linearised equations of 
motion of a rigidly rotating Newtonian polytrope to investi¬ 
gate oscillatory behaviour in the time domain. This will al¬ 
low us to check, as accurately as a time evolution approach 
permits, the f-mode, r-mode and inertial-mode frequencies 
and eigenfunctions in the literature, and also to extend the 
calculation of the inertial mode frequencies into the regime 
of rapid rotation. Our work is in much the same spirit as 
that of Papaloizou & Pringle (1980), who performed time 
evolutions of r-modes in rapidly rotating thin shells. 

Our long-term aim is to provide a stable and well-tested 
platform for developing a more advanced time evolution 
code, capable of evolving more realistic stellar models. For 
instance, future modifications will include evolving a star 
with a differentially rotating background. This issue is of 
particular interest, as recently there has been speculation 
that no mode resembling an r-mode will exist in a sufficiently 
strongly differentially rotating star (Karino et al. 2000). By 
starting with a rigidly rotating configuration, and tracking 
the frequency of a r-mode through a series of runs where 
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the degree of differential rotation is increased step-by-step, 
it should be possible to probe the mode structnre in highly 
differentially rotating stars. 

Another problem that conld be addressed using our 
code is that of modelling gravitational radiation reaction as 
a local force, included in the equations of motion. A number 
of elegant formulations of this problem have been presented 
(see Blanchet et al. 1990 and Rezzolla et al. 1999), but have 
yet to be fully applied to the modes of compact stars. Re¬ 
cently, Lindblom, Tohline & Vallisneri (2001) included radi¬ 
ation reaction in their evolutions of the r-modes using for¬ 
mulae from Rezzolla et al. but made some approximations 
when calculating the force. With our computationally less 
intensive linear code one can attempt to compare the re¬ 
sults obtained via their approximate treatment with those 
obtained from the full formalism, before the mode saturates 
due to nonlinear effects.. 

Another future modification of our code could be to in¬ 
clude the first post-Newtonian corrections to the equilibrium 
star and the equations of motion. This problem is also of con¬ 
siderable interest, as recently Ruoff & Kokkotas (2001) have 
suggested that r-modes will not exist in sufficiently compact 
relativistic neutron stars. They speculate that this is due 
to frame dragging, an effect which is present at first post- 
Newtonian order. This hypothesis could be tested using a 
strategy similar to that described above for differential ro¬ 
tation; specifically by performing a series of evolutions of 
successively more and more compact stars, tracking the r- 
mode between each evolution. 

The implementation of all of these pieces of physics are 
challenging, and we hope to learn how to deal with each 
separately in the linear regime before combining them or 
extending them to the non-linear regime. In essence, our 
strategy is a modular one, with each problem being isolated 
and studied on its own. (As regards non-linear evolutions, 
by coupling several linear time evolution codes it will be 
possible to monitor mode-mode coupling in the weakly non¬ 
linear regime). 

In this paper we postpone such advanced calculations, 
and present evolutions of rigidly rotating Newtonian poly¬ 
tropes, with the perturbations in the gravitational poten¬ 
tial neglected (the Cowling approximation). Our goal is to 
demonstrate the accuracy and stability of our code prior to 
studying the more complicated physical problems discussed 
above. We begin by describing some of the computational 
details (section ^ before examining the f-modes (section 
3.1), th e r- modes (section 3.2) and finally inertial modes 
(section 3.3). 


2 FORMULATION OF THE PROBLEM 
2.1 The unperturbed star 

The unperturbed rotating stellar structure was calculated 
using the Hachisu self-consistent field method (Hachisu 
1986). The equation of state is given by the usual polytropic 
form: 

P^KpJ ( 1 ) 

where P and p denote pressure and density, while K and 7 
are constants. The polytropic exponent 7 is related to the 
polytropic index n by 


7 = 1 - 1 - 1 /n. 


( 2 ) 


Each equilibrium model is characterised by its central den¬ 
sity and ratio of the polar to equatorial axes. In constructing 
an equilibrium model, the equations governing the gravita¬ 
tional potential and the hydrostatic equilibrium are solved 
iteratively in an integral form (starting from a non-rotating 
model) until convergence to a desired numerical accuracy 
is achieved. Our code reproduces the numerical models in 
Hachisu (1986) to the number of significant digits published. 

We will find it useful to introduce a coordinate system 
consisting of the usual spherical polar 0 and (f) coordinates, 
but with a modified radial coordinate x = x{r,9) which is 
fitted to surfaces of constant pressure of the background star. 
With this choice, the unperturbed pressure and density are 
functions of x only, i.e. P = P{x) and p = p{x). The usual 
spherical polar partial derivatives can then be written as: 


d 

dx 

d 

dr 

e dr 

9 dx 

d 

dx 
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( 3 ) 

( 4 ) 


We will choose to set x equal to r in the equatorial plane, 
i.e. a;(r, 7 r/ 2 ) = r. For 9 < 7 r/ 2 , the oblate shape of the star 
will mean that x(r,9) > r. The partial derivatives of x with 
respect to r and 9 can be easily evaluated using: 
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e dr 
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( 5 ) 

( 6 ) 


This choice of coordinate system will make it particularly 
easy to im plem ent the outer boundary conditions described 
in section (2.3). 


2.2 Equations to be solved 

For a given background star, three pieces of physics are re¬ 
quired to describe this problem: The equations of motion, 
the equation of mass conservation, and an equation of state 
for the perturbations. 

The equations of motion referred to the frame of the 
rotating background are: 

dv 

— -I- (v • V)v -I- 2n X V -I- n X (n X r) 

= _ivp-v<(., (7) 

p 

where v denotes the departure of the fluid from the rigid 
rotation of the background equilibrium model. Performing 
an Eulerian perturbation and making the Cowling approxi¬ 
mation leads to: 

^ -f 2H X V = %VP - -V5P. (8) 

ot p 

The prefix 5 denotes an Eulerian perturbation, while p and 
P refer to the background star. The equation of continuity 
is: 

|^+V.(pv)=0 (9) 

which linearises to give: 
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f: + v.(,v) = o. 


( 10 ) 


We will write the equation of state for the perturbations as: 


AP _ Ap 

—p- — 7pert , 


( 11 ) 


where A denotes a Lagrangian perturbation. Using the re¬ 
lation that, for a displacement vector the Lagrangian and 
Eulerian perturbations of a scalar are related by 


A = 5-bC- V 

(12) 

we find: 


' ‘‘f e-A. 

(13) 

P 7pert P 

The vector A is given by: 


A = V log p-^—V log P, 

(14) 



and its magnitude is known as the Schwarzschild discrimi¬ 
nant, which depends upon the background pressure and den¬ 
sity, and on the index 7pert describing the perturbations. We 
will consider the case where the perturbations obey the same 
equation of state as the background, so that 7pert = 7 and 


A = 0. 


(15) 


As discussed in the stellar oscillation literature (see e.g. Tas- 
soul [1978] or Unno et al. [1989]), the magnitude and direc¬ 
tion of A determines the nature of g-modes in the star. If A 
points inwards the star has stable oscillatory g-modes, while 
if A points outwards the g-modes are unstable, and the star 
is described as being ‘convectively unstable’. In choosing 
A = 0 there exists no restoring force at all for g-mode type 
displacements, so we have eliminated the g-modes altogether 
from the stellar spectrum. 

Expressed in terms of SP, equations (^) and (p^ be- 


come: 




VP ISP 1 r 

(16) 

—h 2S7 XV = 
at 

- --VSP, 

p 7 P p 

dSP 

dt 

-^V.(pv). 

(17) 


With respect to the standard spherical polar vector basis 
[er,ee,e,^] {not the coordinate vector basis) we then find: 


dfr 

dt 


-f 2flsinf 
or 


'> "o- d" 

or j P 


dfe 

dt 


13(5P IdPlSP 

- +2ilcos9f^ + 

r oO r oO j P 


(18) 

(19) 


df± 

dt 


1 dSP 
r sin 6 dcp 


2fI(cos dfg + sin dfr), 


( 20 ) 


dSP 

dt 



f. 


( 21 ) 


where the velocity v has been replaced by the flux f = pv, as 
this choice of variable lead s to very simple outer boundary 
conditions (see section [ 2 !^ . 

We follow Papaloizou & Pringle (1980) in decomposing 
any given perturbation as a sum over basis functions of the 
form cosm(f) and sinmi^i, where m is an integer, e.g. for the 
pressure: 


m = + oo 

SP{t,r,e,(f)) = E SPf {t, r, 9) cos m(f) 

m=0 


-b 5Pm {t, r, 9) sin mij). (22) 

This allows us to reduce the number of space dimen¬ 
sions that appear in the numerical evolution from three to 
two. With the understanding that the perturbation vari¬ 
ables are now functions of {t,r,9), but not rj), equations 
(@)-(il) then apply to each azimuthal number m sep¬ 
arately, i.e. the equations decouple in m, giving a total 
of eight equations. The partial derivatives with respect to 
(f are trivial to calculate analytically. In summary, equa¬ 
tions (p^)-(pl|) then give a set of eight coupled first order 
linear partial differential equations in the eight unknowns 
{5P+,5P-, /+ , fr- . /+ , U . ff Jf)- 

For a perturbation of definite azimuthal number m, this 
decomposition with respect to (f> greatly reduces the mem¬ 
ory storage and running time of the code, as compared to a 
fully three-dimensional computation. The fundamental rea¬ 
son that the equations decoupled is that the background star 
is axisymmetric. A further decomposition in which the 9 de¬ 
pendence of the perturbation variables is written in terms of 
a set of spherical harmonic functions Yim is not very useful— 
the non-spherical nature of the rotating unperturbed star 
would couple terms with different I indices, with the cou¬ 
pling becoming stronger the more rapidly the star rotates. 
By eliminating only the 4> dependence, we have simplified 
the computational costs of the problem to the maximum ex¬ 
tend possible without significantly increasing the complexity 
of the equations to be solved. 


2.3 Boundary conditions 

By definition, the Lagrangian perturbation in the pressure 
is zero at the surface, i.e. 

SP + ^ ■ \/P = 0 (at surface). 

For the unperturbed star, equation (Q) gives: 

VP = —pfl X (n X r) — p\7(j). (23) 

For a polytrope p = 0 at the surface, from which it follows 
that VP = 0 at the surface. It follows at once that our outer 
boundary condition is: 

SP = 0 (at surface). (24) 

Our choice of coordinates makes imposing the outer bound¬ 
ary condition very simple: We have SP{x = R) — 0, where 
R denotes the equatorial radius. Also, the mass flux f = pv 
is zero at the surface, by virtue of the vanishing density, and 
so 

f = 0 (at surface). (25) 

Our choice of surface fitted coordinates and perturbation 
variables guarantees that we never need to evolve quantities 
at the surface of the star. Numerically, this is highly desir¬ 
able, as resolving the possibly rapidly falling density profile 
near the surface can lead to inaccuracies and instabilities. 

As we will consider perturbations with m > 0 we can 
set f = 0 at the centre. 


3 RESULTS 

Equations (^^-(^^, dec omp osed with respect to (j) as de¬ 
scribed above (equation fc^) were propagated forward in 
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Figure 1. Demonstration of second order convergence. The three 
upper curves plot the kinetic energy at three different grid reso¬ 
lutions, labelled as high, medium and low resolutions. The lower 
curve is the dimensionless convergence ratio, defined in equation 
(bm. Time is plotted in unity of l/y'TrGpoi where po is the average 
density of the non-rotating star of the same mass. 


time using the two-step Lax-WendrofF scheme described in 
Press et al. (1986). The code is second order convergent, 
and runs stably for hundreds of rotation periods. The sec¬ 
ond order convergence is illustrated by figure In this fig¬ 
ure we plot the total kinetic energy of the star, as measured 
in the rotating frame, as a function of time for some arbi¬ 
trary initial data. The time is given in units of llJJGpo, 
where po is the average density of the non-rotating star of the 
same mass and equation of state. The star rotates at a rate 
Q/y^nGpo = 0.538, corresponding to a polar-to-equatorial 
axis ratio of 0.79. The evolution was performed for three 
different grid resolutions, labelled as high, medium and low 
in the figure, with the resolution differing by a factor of two 
between each successive run. In full, for the high resolution 
run we used a grid of 64 angular points by 60 radial ones, 
the medium resolution grid was 32 by 30 points, and the 
low resolution grid was 16 by 15 points. The kinetic energy 
is normalised to unity at t = 0. A small oscillation is visi¬ 
ble, corresponding to the excitation of a p-mode (Cox 1980), 
whose exact nature need not concern us here. A convergence 
ratio C is then calculated as the ratio: 

i?(t, high resolution) — i?(t, low resolution) 
if (t, high resolution) — i?(t, medium resolution) 

For a code with second order accuracy in space, the ratio 
can be shown to be approximately equal to 5. This is indeed 
the case for our time evolution code. Beyond the time inter¬ 
val shown in the figure, the oscillations in the energy start 
to loose phase coherence, and the convergence ratio begins 
to oscillate strongly. The convergence ratio then looses its 
significance. 

The long-term stability of the code is illustrated by fig¬ 
ure]^, which plots the kinetic energy of the star as a function 
of time. The rotation rate is the same as above, and the star 
rotates 100 times over the period plotted. A small amount of 
artificial viscosity was included to give this long-term stabil¬ 


Figure 2. Kinetic energy as a function of time, normalised to 
nity at t = 0, and time measured in the same units as in figure 
). Around 100 stellar rotations occur over the period plotted. 

ity, and the grid consisted of 32 angular points and 30 radial 
ones. The initial data used was that of an Z = m = 2 r-mode 
for a slowly rotating star (about 80 r-mode oscillations oc¬ 
cur over this interval). A rather sharp loss of kinetic energy 
occurs over the first ten or so rotations. This is likely due to 
the initial data exciting a number of short-wavelength modes 
whose energy was rapidly dissipated. Over the remainder of 
the evolution, the kinetic energy remains approximately con¬ 
stant, falling by about 10% over the remaining 90 rotations, 
confirming the long-term stability of our code. This evolu¬ 
tion required approximately two hours of computer time on 
a standard PC. Most of the results in this paper are based 
on runs of this duration or less. 

Part of the above energy loss is due to the artificial vis¬ 
cosity, and part due to numerical error. Given that the code 
has been shown to be second order convergent, an increase 
in grid resolution by a factor of N would reduce the rate 
of energy loss due to numerical error by a factor On 
the basis of runs we have performed where the artifical vis¬ 
cosity was set to zero, we estimate that an increase in grid 
resolution by a factor of 10 (i.e. a 320 by 300 grid) would 
be sufficient to reduce the energy loss to around 1% for an 
evolution of the above duration. Such a run would take over 
a week on a standard PC, and so has not been performed 
here. However, should such a high level of accuracy be re¬ 
quired for a particular application, the necessary runs could 
be easily be done on a supercomputer. 

In the following sections we will use this time evolution 
code to study f-modes, r-modes and inertial modes. We will 
use as initial data some (possible crude) approximation to 
the true eigenfunction. The perturbed quantities will then 
be Fourier analysed. By performing a series of such evolu¬ 
tions, of gradually increasing stellar rotation rate, the mode 
of interest can then be tracked from the slow-rotation case 
(where it’s approximate frequency is already known) into 
the rapid rotation regime. 

Note that for an evolution of duration T, the error in 
the frequency localisation of the transform is approximately 
equal to 1/T. It was therefore necessary to perform long 
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time evolutions to obtain accurate frequencies. This error 
in measurement is indicated by error bars on the mode fre¬ 
quency plots that follow (figures However, it should 

be noted that this is not the only source of error—the finite 
grid resolution will affect the results also, but its effect is 
not as easy to quantify. The error bars in the figures should 
be regarded as simple estimates of the error. 

Note that all mode frequencies given in this paper are 
measured in the rotating frame. Also, we only present results 
for the case when the polytropic index n is equal to unity. 

3.1 The f-modes 

As described above, to investigate a particular mode it was 
necessary to supply our code with initial data that would 
excite it. In the case of f-modes we chose to use initial data 
of the form: 

where R{0) denotes the stellar radius at polar angle 9 and 
Yu is a spherical harmonic. This reduces to the correct f- 
mode pressure perturbation in the limit of zero rotation for 
an incompressible star, in the Cowling approximation. The 
perturbed mass flux f was set to zero in the initial data. 
The mode frequency was then extracted by taking power 
spectra of the time series of the perturbed quantities. A 
more efficient excitation of the f-modes would have been 
possible if non-zero initial data was provided for f also, but 
this was not found to be necessary—the f-modes were easily 
identifiable using the above initial data. Indeed, this scheme 
efficiently excited the I = m mode in all but the most rapidly 
rotating stars. 

The power spectra of the evolutions consisted of just one 
peak for the non-rotating star. This then split into two peaks 
for 7 ^^ 0 , corresponding to the breaking of the degeneracy 
between co- and counter-rotating modes. This is illustrated 
in figure ^ for the spectra of SP'^ for the Z = m = 2 f-mode. 
The left hand panel was obtained for a non-rotating star, 
while the right hand one was for a rotating one (fl/y^wGpo = 
0.37, polar-to-equatorial radius ratio 0.92). Note how the 
power in the degenerate = 0 peak for the static star is 
split roughly equally between the co- and counter-rotating 
peaks in the 7 ^ 0 case. As the star’s rotation rate was 
increased, the frequency separation of the modes increased, 
and the power was shared less evenly between the two. 

Before examining the case of rotating stars in detail 
we will consider the simpler case of non-rotating ones. The 
mode frequencies, as calculated by a normal mode analysis 
in the Cowling approximatio n (K okkotas, personal commu¬ 
nication), are given in table |3.l| for the / = m = 2, 4 and 
6 modes. Also shown are the mode frequencies we obtained 
from our time evolution code, and the percentage difference 
between the two. In this non-rotating case the eigenmode 
calculation is known to be highly accurate, so that the per¬ 
centage differences in the table give an indication of the 
accuracy of our code. The difference is very small (less than 
1 %) in the I = m = 2 case, but somewhat higher (around 
4%) for the high-index modes. This decrease in accuracy 
is almost certainly due the higher-index modes varying on 
smaller length scales than the Z = m = 2 mode, and there¬ 
fore being less well represented on our spatial grid. 


Non-rotating star Rotating star 



Figure 3. Power spectra of the perturbed pressure (5P+ for a non¬ 
rotating star (left panel) and for a rotating star (right panel). The 
powers are in the same (arbitrary) units in both panels, while the 
mode frequencies are in units of y/irGpo, as indicated. 


1 = m = 

2 

4 

6 

Eigenmode 

1.912 

2.577 

3.077 

Time evolution 

1.90 

2.68 

2.96 

Percentage difference 

0.5 

3.8 

3.9 


Table 1. F-mode frequencies as calculated by an eigenmode anal- 
sis, and using our time evolution code (units same as in figure 
). The eigenmode vales are anticipated to be of high precision, 
so that the percentage difference between the two is an indication 
of the accuracy of our time evolutions. 


We then investigated f-modes in rotating stars. The f- 
mode eigenfrequencies and eigenfunctions of rapidly rotat¬ 
ing stars have been calculated by Ipser & Lindblom (1990), 
without making the Cowling approximation. Specifically, 
they considered polytropic stars, and searched for modes 
which had spherical harmonic indices I = m in the limit of 
slow rotation. 

The frequency of the co- and counter-rotating modes, 
together with the Ipser & Lindblom values for the counter¬ 
rotating mode, are plotted in figures ^ and for I = m = 
2, 4, 6 respectively. 

For most rotation rates the counter-rotating mode fre¬ 
quencies calculated from our time evolutions are somewhat 
higher than those calculated by Ipser & Lindblom. This is al¬ 
most certainly due to our having made the Cowling approx¬ 
imation, while Ipser & Lindblom have not. The discrepancy 
is typically of order 30% for I = m = 2, 10% for Z = m = 4 
and 3% for Z = m = 6 . Such a decrease of discrepancy with 
increasing harmonic index is to be expected: higher order 
modes have alternating regions of positive and negative den¬ 
sity perturbation, whose contribution to the perturbation in 
the gravitational potential at a given point tend to cancel 
one-another out (Cowling 1941). We can therefore conclude 
that, within the error margin of the Cowling approximation, 
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Figure 4. The I = m = 2 f-mode frequencies as a function of 
rotation rate. The upper curve is for the counter-rotating mode 
as calculated by our time evolution code. This is to be compared 
with the Ipser & Lindblom result of the solid curve. The dif¬ 
ference between the two results is due to our having made the 
Cowling approximation. The monotonically decreasing curve is 
the co-rotating mode frequency as calculated by our code. All 
frequencies are in units of y/nGpo, where po is the average den¬ 
sity of the non-rotating star of the same mass and equation of 
state. In these units, the Kepler rotation rate is 2/3. 



iJ/V(ltGpj) 

Figure 5. As figure ^ with I = m = A 

our time evolution code confirms the f-mode frequencies cal¬ 
culated by Ipser & Lindblom. 


3.2 The r-modes 

Because of the considerable interest aroused by the r-modes 
(see Andersson & Kokkotas 2001 for a recent review), we 
used our time evolution code to examine both the frequen¬ 
cies and eigenfunctions of this class of mode. For initial data 
we provided a perturbation in the mass flux of the form: 



Figure 6. As figure ^ with I = m = 6 



where Y^n is a magnetic spherical harmonic (Thorne 1980). 
In the limit of slowly rotating stars, this becomes an exact 
mode solution of the linearised equations of motion. 


3.2.1 The r-mode frequencies 

In order to extract r-mode frequencies, the power spectra 
of the 9 and (j> components of the mass flux were analysed. 
The result of a series of such evolutions for the I = m = 2 
r=mode is shown in figure |^, with the error bars indicating 
the uncertainty in the frequency. On the vertical axis we plot 
the dimensionless ratio cr/fl, where a is the mode frequency. 
On the horizontal axis we plot the normalised rotation rate 
Q,/y/JGpo. For comparison, a curve showing the mode fre¬ 
quency as calculated by Karino et al. (2000) is shown also, 
together with a horizontal line at the mode frequency as 
calculated using the slow rotation limit {a/Q. = 2/3). Note 
that Karino et al. did not make the Cowling approximation. 
Within the indicated uncertainties due to the finite duration 
of the time evolution, the data sets agree very well, even at 
the fastest rotation rates. To leading order in the rotation 
rate, the Cowling approximation is known to not affect the 
mode frequency, so this close agreement was to be expected. 
There is a slight difference between our results and those 
of Karino et al. at the highest rotation rates. We suspect 
this is due to the Cowling approximation being a poorer 
approximation in this regime, as the fractional density per¬ 
turbation 5p/p is known to increase with rotation rate as 
for r-modes (Andersson & Kokkotas 2001). 


3.2.2 The r-mode eigenfunctions 

The r-mode eigenfunctions were calculated by Karino et al. 
(2000). In order to confirm their veracity, these eigenfunc¬ 
tions were used as initial data for our time evolutions. An 
accurate eigenfunction would excite the appropriate r-mode 
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— Slow rotalion limit (analytic) 

— Karino et al. 

Time evolution code 
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Figure 7. The frequency of the I = m = 2 r-mode as a function of 
rotation rate. Mode frequencies are given in units of Q, the stellar 
rotation rate. As in previous plots, the error bars indicated the 
result of our time evolution code. The dotted line gives the mode 
frequency calculated by Karino et al. (2000). The horizontal line 
indicates the analytic slow rotation limit of 2/3. 


(t/ ^nGpo = 

0.170 

0.538 

0.624 



Peak ratio 


N 

224 

42.6 

4.26 

5P+ 

7740 

174 

8.06 


Table 2. Ratio of m = 2 r-mode excitation peak to next strongest 
peak (representing other modes), for the quantities N <SP+, 

for three different rotation rates. Recall than the Kepler rotation 
rate is 2/3 in these units. 


only, with little excitation of other modes, while an inaccu¬ 
rate function would deposit a signihcant quantity of energy 
over other modes. In order to quantify this, we examined 
the power spectra of several evolutions with different stellar 
rotation rates, for the I = m = 2 r-mode. For each evolution 
we noted the ratio of the r-mode excitation peak to the next 
strongest excitation. 

We found that in the power spectra of fe and f^, the 
r-mode excitation was extremely clean, so that secondary 
peaks are many orders of magnitude smaller. For fr and Sp 
significant secondary peaks were present, implying that our 
initial data also excited some radial modes or polar modes 
of oscillation. The ratios of these to the r-mode peak are 
given in table 3.2.2 , for stars with rotation rates Q/y^nGpo 
of 0.170, 0.538 and 0.624. 


Clearly, in all cases the initial data efficiently picks out 
the r-mode from the stellar spectrum. The excitation is more 
accurate the slower the rotation of the star, but even in the 
the rapidly rotating fl/y^TrUpo = 0.624 case there is little 
other mode excitation. The kinetic energy of the perturba¬ 
tions is quadratic in the mass flux, so that in this case only 
a few percent of the kinetic energy of the initial data is de¬ 
posited into modes other than the r-mode. For the relatively 


slowly rotating ^/^nGpo = 0.170 case, only a few thou¬ 
sandths of a percent is deposited. We can therefore confirm 
that the eigenfunctions computed by Karino et al. (2000) 
are accurate. 


3.3 Inertial modes 

The r-modes of the previous section are a special case 
of a wider class of oscillations, known as inertial modes 
(Greenspan 1968, Pedlosky 1987). The dehning feature of 
the inertial modes is that their frequency vanishes linearly 
in the limit of zero stellar rotation, or equivalently that cr/fl 
remains Hnite in this limit. The r-modes of the last section 
are the subset of the inertial modes which have purely ax¬ 
ial velocity perturbations in the limit of zero rotation (axial 
means that the velocity can be written in terms of V x Tim). 
Furthermore, for stars with zero Schwarzchild discriminant, 
the two spherical harmonic indices of the r-modes must be 
equal, i.e. I = m (Provost et al. 1981). 

In contrast, inertial modes are hybrid, having veloc¬ 
ity perturbations with both axial and polar parts, even in 
the slow rotation limit (polar means that the velocity field 
can be written as a sum of terms proportional to YirnSr 
and VYim)- Also, the inertial modes are not restricted to 
I = m spherical harmonics. Inertial modes were considered 
in detail recently by Lockitch & Friedman (1999), who cal¬ 
culated mode frequencies and eigenfunctions. It was found 
that many of these inertial modes could undergo the CFS 
instability and are therefore interesting from a gravitational 
wave point of view, even though they are probably not as 
strongly unstable as the I = m = 2 r-mode. It is therefore 
of interest to extend the calculations of Lockitch & Fried¬ 
man and calculate these mode frequencies for more rapidly 
rotating stars. 

We have performed such an analysis using our time evo¬ 
lution code. We used the analytic approximations to eigen¬ 
functions given by Lockitch & Friedman as initial data to 
carry out evolutions of several inertial modes. 

A graph showing the frequency of three such modes as 
a function of stellar rotation rate is shown in hgure These 
are m = 2 modes whose leading contribution to the velocity 
held has harmonic index I — 4. Analytic hts to the radial 
eigenfunctions of these modes in uniform density stars are 
given in table 3 of Lockitch & Friedman. As the latter au¬ 
thors note, the eigenfunctions are very similar for n = 1 
polytropes, so we have used these analytic functions to sup¬ 
ply initial data to our code. As in previous hgures, the verti¬ 
cal axis plots the dimensionless mode frequency as measured 
in the rotating frame, while the angular rotation rate is given 
in units of y/nGpo. The three horizontal lines indicate the 
mode frequency in the slow rotation limit as calculated by 
Lockitch & Friedman (see table 6 of their paper). In the limit 
of slow rotation we found the eigenfunctions were sufficiently 
accurate to deposit approximately 99% of their energy into 
a single mode, with frequency very close to that computed 
by Lockitch & Friedman, conhrming our identihcation of the 
mode. The slight difference between the frequencies is due to 
our having made the Cowling approximation, while Lockitch 
& Friedman have not. 
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Figure 8. Oscillation frequencies for inertial modes. The three 
modes all have leading order velocity perturbations proportional 
to the / = 4, m = 2 spherical harmonic in the limit of slow rota¬ 
tion, but have different radial behaviour. The units are the same 
as in figure M. The horizontal lines indicate the mode frequen¬ 
cies as calculated by Lockitch & Friedman in the slow rotation 
non-Cowling limit. 


4 CONCLUSIONS 

In this paper we have presented results from a numerical 
code for time evolutions of perturbations of rapidly and 
rigidly rotating Newtonian polytropes (within the Cowling 
approximation). Our code is second order convergent and 
runs stably for many (i.e. hundreds) of stellar oscillations. 
We have compared the results of our evolutions with calcula¬ 
tions of f- and r-modes in the literature and found agreement 
within the errors expected for the Cowling approximation. 
We have also presented results for the inertial modes, and 
(for the first time) extended the frequencies calculated by 
Lockitch & Friedman into the rapid rotation regime. 

Even though we have provided some new results for 
oscillations in stars rotating near the break-up limit, this 
work is essentially a “proof of principle”. Our main aim was 
to illustrate the accuracy attainable with a time-evolution 
approach to the problem of rotating stars. In particular, the 
relative simplicity and reliability of the perturbative 2D ap¬ 
proach makes it a valuable complement to fully 3D hydrody- 
namical simulations, such as the recent ones by Lindblom, 
Tohline and Vallisneri (2001). 

Having demonstrated the stability and accuracy of our 
code, we plan to use it as a test-bed for the implementation 
of additional physical features, which we will add in a mod¬ 
ular fashion, one at a time. Via such a gradual escalation in 
complexity, we hope to arrive at a code capable of evolving 
rather more realistic stars, incorporating (for instance) dif¬ 
ferential rotation, gravitational radiation reaction, and non¬ 
linear effects. Such extensions are currently underway, and 
will be presented in due course. 
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